#pragma once
#include "../Utility/Log.hpp"
#include "Vector2.hpp"
#include "Vector.hpp"
#include "../common.hpp"

namespace zzz{
template <unsigned int ROW, unsigned int COL, typename T>
class MatrixBase
{
public:
  // type definitions
  typedef T               value_type;
  typedef T*              iterator;
  typedef const T*        const_iterator;
  typedef T&              reference;
  typedef const T&        const_reference;
  typedef zuint           size_type;
  typedef size_type       difference_type;

protected:
  Vector<ROW,Vector<COL, T> > V;

public:
  //STL support
  // iterator support
  iterator begin() { return (T*)&V; }
  const_iterator begin() const { return Data(); }
  iterator end() { return Data()+ROW*COL; }
  const_iterator end() const { return Data()+ROW*COL; }

  // reverse iterator support
  typedef std::reverse_iterator<iterator> reverse_iterator;
  typedef std::reverse_iterator<const_iterator> const_reverse_iterator;

  reverse_iterator rbegin() { return reverse_iterator(end()); }
  const_reverse_iterator rbegin() const { return const_reverse_iterator(end());}
  reverse_iterator rend() { return reverse_iterator(begin()); }
  const_reverse_iterator rend() const {return const_reverse_iterator(begin());}

  // at() with range check
  const T& at(const size_type pos) const {
    ZRCHECK_LT(pos, ROW*COL);
    return Data()[pos];
  }
  T& at(const size_type pos) {
    ZRCHECK_LT(pos, ROW*COL);
    return Data()[pos];
  }

  // front() and back()
  reference front() { return V[0][0]; }
  const_reference front() const {  return V[0][0];}
  reference back() {   return V[ROW][COL]; }
  const_reference back() const {   return V[ROW][COL]; }

  // size is constant
  static size_type size() { return ROW*COL; }
  static bool empty() { return false; }
  static size_type max_size() { return ROW*COL; }
  enum { static_size=ROW*COL };

  // swap (note: linear complexity)
  void swap (MatrixBase<ROW,COL, T>& y) {  std::swap_ranges(begin(),end(), y.begin());}

  //subscript
  const T& operator[](const size_type pos) const {return at(pos);}
  T& operator[](const size_type pos) {return at(pos);}
  
  const T& At(const size_type x, const size_type y) const {
    ZRCHECK_LT(x, ROW);
    ZRCHECK_LT(y, COL);
    return V[x][y];
  }
  T& At(const size_type x, const size_type y) {
    ZRCHECK_LT(x, ROW);
    ZRCHECK_LT(y, COL);
    return V[x][y];
  }
  const T& At(const Vector<2,size_type> &pos) const {return At(pos[0],pos[1]);}
  T& At(const Vector<2,size_type> &pos) {return At(pos[0],pos[1]);}
  
  const T& operator()(const size_type x, const size_type y) const {return At(x, y);}
  T& operator()(const size_type x, const size_type y) {return At(x, y);}
  const T& operator()(const Vector<2,size_type> &pos) const {return At(pos[0],pos[1]);}
  T& operator()(const Vector<2,size_type> &pos) {return At(pos[0],pos[1]);}

  const Vector<COL, T>& Row(const size_type row) const{
    ZRCHECK_LT(row, ROW);
    return V[row];
  }
  Vector<COL, T>& Row(const size_type row) {
    ZRCHECK_LT(row, ROW);
    return V[row];
  }

  //copy
  const MatrixBase<ROW,COL, T> &operator=(const T &a){for (size_type i=0; i<ROW*COL; i++) Data()[i]=a; return *this;}
  const MatrixBase<ROW,COL, T> &operator=(const MatrixBase<ROW,COL, T> &a){IOObject<T>::CopyData(Data(), a.Data(),ROW*COL); return *this;}
  template<typename T1>
  const MatrixBase<ROW,COL, T> &operator=(const MatrixBase<ROW,COL,T1> &a){for (zuint i=0; i<ROW*COL; i++) Data()[i]=a[i]; return *this;}
  
  //raw data access
  inline T* Data(){return reinterpret_cast<T*>(&V);}
  inline const T* Data() const {return reinterpret_cast<const T*>(&V);}
  
  //rawcopy
  void RawCopy(const MatrixBase<ROW,COL, T> &other){memcpy(Data(),other.Data(),sizeof(T)*ROW*COL);}
  void RawCopy(const T *data){memcpy(Data(),data,sizeof(T)*ROW*COL);}
  void Zero(zuchar x=0) {memset(Data(), x, sizeof(T)*ROW*COL);}
  static MatrixBase<ROW,COL, T> GetZero(int x=0) {MatrixBase<ROW,COL, T> mat; mat.Zero(x); return mat;}
  void Fill(const T &a){for (size_type i=0; i<ROW*COL; i++) Data()[i]=a;}
  static MatrixBase<ROW,COL, T> GetZero(const T x) {MatrixBase<ROW,COL, T> mat; mat.Fill(x); return mat;}

  //////////////////////////////////////////////////////////////////////////
  // index
  size_type ToIndex(const Vector<2,size_type> &pos) const {return pos[0]*COL+pos[1];}
  Vector<2,size_type> ToIndex(size_type idx) const {
    Vector<2,size_type> pos;
    pos[0]=idx/COL;pos[1]=idx-pos[0]*COL;
    return pos;
  }
  bool CheckIndex(const Vector<2,size_type> &pos) const {return (pos[0]<ROW && pos[1]<COL);}
  bool CheckIndex(const size_type &pos) const {return pos<ROW*COL;}

  //Min and Max
  inline T Max() {
    T maxv=Data()[0];
    for (unsigned int i=1; i<ROW*COL; i++) if (maxv<Data()[i]) maxv=Data()[i];
    return maxv;
  }

  inline int MaxPos() {
    T maxv=Data()[0];
    int pos=0;
    for (unsigned int i=1; i<ROW*COL; i++) if (maxv<Data()[i]) {maxv=Data()[i];pos=i;}
    return pos;
  }

  inline T Min() {
    T minv=Data()[0];
    for (unsigned int i=1; i<ROW*COL; i++) if (minv>Data()[i]) minv=Data()[i];
    return minv;
  }

  inline int MinPos() {
    T minv=Data()[0];
    int pos=0;
    for (unsigned int i=1; i<ROW*COL; i++) if (minv>Data()[i]) {minv=Data()[i];pos=i;}
    return pos;
  }

  inline T AbsMax() {
    T maxv=Data()[0];
    for (unsigned int i=1; i<ROW*COL; i++) if (maxv<abs(Data()[i])) maxv=Data()[i];
    return maxv;
  }

  inline int AbsMaxPos() {
    T maxv=Data()[0];
    int pos=0;
    for (unsigned int i=1; i<ROW*COL; i++) if (maxv<abs(Data()[i])) {maxv=Data()[i];pos=i;}
    return pos;
  }

  inline T AbsMin() {
    T minv=Data()[0];
    for (unsigned int i=1; i<ROW*COL; i++) if (minv>abs(Data()[i])) minv=Data()[i];
    return minv;
  }

  inline int AbsMinPos() {
    T minv=Data()[0];
    int pos=0;
    for (unsigned int i=1; i<ROW*COL; i++) if (minv>abs(Data()[i])) {minv=Data()[i];pos=i;}
    return pos;
  }

  // array math
  inline T Dot(const MatrixBase<ROW,COL, T> &a) const {T ret=0;for (unsigned int i=0; i<ROW*COL; i++) ret+=Data()[i]*a.Data()[i]; return ret;}
  inline T Sum() const{T ret=0;for (unsigned int i=0; i<ROW*COL; i++)  ret+=Data()[i];return ret;}
  inline void Negative() {for (size_type i=0; i<ROW*COL; i++) Data()[i]=-Data()[i];}
  inline T Len() const {return (T)sqrt((double)LenSqr());}
  inline T LenSqr() const {return Dot(*this);}

  inline T Normalize() {T norm=Len(); *this/=norm; return norm;}
  inline MatrixBase<ROW,COL, T> Normalized() const {return *this/Len();}

  inline T DistTo(const MatrixBase<ROW,COL, T> &a) const {return (T)sqrt((double)DistToSqr(a));}
  inline T DistToSqr(const MatrixBase<ROW,COL, T> &a) const {MatrixBase<ROW,COL, T> diff(*this-a); return diff.Dot(diff);}

  // comparisons
  bool operator== (const MatrixBase<ROW,COL, T>& y) const {return std::equal(begin(), end(), y.begin());}
  bool operator< (const MatrixBase<ROW,COL, T>& y) const {return std::lexicographical_compare(begin(),end(), y.begin(), y.end());}
  bool operator!= (const MatrixBase<ROW,COL, T>& y) const {return !(*this==y);}
  bool operator> (const MatrixBase<ROW,COL, T>& y) const {return y<*this;}
  bool operator<= (const MatrixBase<ROW,COL, T>& y) const {return !(y<*this);}
  bool operator>= (const MatrixBase<ROW,COL, T>& y) const {return !(*this<y);}  

  inline bool operator==(const T &c) const {
    for (int i=0; i<ROW*COL; i++)
      if (Data()[i]!=c) return false;
    return true;
  }

  inline bool operator!=(const T &c) const {return !(*this==c);}

  //IO
  friend inline ostream& operator<<(ostream& os,const MatrixBase<ROW,COL, T> &me) {
    for (int i=0; i<ROW; i++) {
      for (int j=0; j<COL; j++)
        os<<me.V[i][j]<<' ';
      if (i!=ROW-1) os<<'\n';
    }
    return os;
  }

  friend inline istream& operator>>(istream& is,MatrixBase<ROW,COL, T> &me) {
    for (int i=0; i<ROW*COL; i++)
      is>>me.Data()[i];
    return is;
  }

  void SaveToFileA(ostream &fo) {
    for (int i=0; i<ROW*COL; i++)
      fo<<Data()[i]<<' ';
  }

  void LoadFromFileA(istream &fi) {
    for (int i=0; i<ROW*COL; i++)
      fi>>Data()[i];
  }

  //positive and negative
  inline const MatrixBase<ROW,COL, T> operator+() const {return *this;}
  inline const MatrixBase<ROW,COL, T> operator-() const {
    MatrixBase<ROW,COL, T> ret;
    for (int i=0; i<ROW*COL; i++)
      ret.Data()[i]=-Data()[i];
    return ret;
  }

  //addition
  inline MatrixBase<ROW,COL, T> operator+(const MatrixBase<ROW,COL, T> &a) const {
    MatrixBase<ROW,COL, T> ret;
    for (unsigned int i=0; i<ROW*COL; i++)
      ret[i]=Data()[i]+a.Data()[i];
    return ret;
  }

  inline MatrixBase<ROW,COL, T> operator+(const T &a) const {
    MatrixBase<ROW,COL, T> ret;
    for (unsigned int i=0; i<ROW*COL; i++)
      ret[i]=Data()[i]+a;
    return ret;
  }

  inline void operator+=(const MatrixBase<ROW,COL, T> &a) {
    for (unsigned int i=0; i<ROW*COL; i++)
      Data()[i]+=a[i];
  }

  inline void operator+=(const T &a) {
    for (unsigned int i=0; i<ROW*COL; i++)
      Data()[i]+=a;
  }

  //subtraction
  inline MatrixBase<ROW,COL, T> operator-(const MatrixBase<ROW,COL, T> &a) const {
    MatrixBase<ROW,COL, T> ret;
    for (unsigned int i=0; i<ROW*COL; i++)
      ret.Data()[i]=Data()[i]-a.Data()[i];
    return ret;
  }

  inline MatrixBase<ROW,COL, T> operator-(const T &a) const {
    MatrixBase<ROW,COL, T> ret;
    for (unsigned int i=0; i<ROW*COL; i++)
      ret.Data()[i]=Data()[i]-a;
    return ret;
  }

  inline void operator-=(const MatrixBase<ROW,COL, T> &a) {
    for (unsigned int i=0; i<ROW*COL; i++)
      Data()[i]-=a.Data()[i];
  }

  inline void operator-=(const T &a) {
    for (unsigned int i=0; i<ROW*COL; i++)
      Data()[i]-=a;
  }

  //matrix multiplication
  template<unsigned int COL2>
  inline MatrixBase<ROW,COL2, T> operator*(const MatrixBase<COL,COL2, T> &a) const {
    MatrixBase<ROW,COL2, T> ret;
    ret.Zero();
    for (int i=0; i<ROW; i++) for (int j=0; j<COL2; j++) {
      for (int k=0; k<COL; k++)
        ret.At(i, j)+=At(i, k)*a.At(k, j);
    }
    return ret;
  }

  //multiply
  inline MatrixBase<ROW,COL, T> ScaleMultiply(const MatrixBase<ROW,COL, T> &a) const {
    MatrixBase<ROW,COL, T> ret;
    for (unsigned int i=0; i<ROW*COL; i++)
      ret.Data()[i]=Data()[i]*a.Data()[i];
    return ret;
  }

  inline MatrixBase<ROW,COL, T> operator*(const T &a) const {
    MatrixBase<ROW,COL, T> ret;
    for (unsigned int i=0; i<ROW*COL; i++)
      ret.Data()[i]=Data()[i]*a;
    return ret;
  }

  inline void operator*=(const MatrixBase<ROW,COL, T> &a) {
    for (unsigned int i=0; i<ROW*COL; i++)
      Data()[i]*=a.Data()[i];
  }

  inline void operator*=(const T &a) {
    for (unsigned int i=0; i<ROW*COL; i++)
      Data()[i]*=a;
  }

  //division
  inline MatrixBase<ROW,COL, T> operator/(const MatrixBase<ROW,COL, T> &a) const {
    MatrixBase<ROW,COL, T> ret;
    for (unsigned int i=0; i<ROW*COL; i++)
      ret.Data()[i]=Data()[i]/a.Data()[i];
    return ret;
  }

  inline MatrixBase<ROW,COL, T> operator/(const T &a) const {
    MatrixBase<ROW,COL, T> ret;
    for (unsigned int i=0; i<ROW*COL; i++)
      ret.Data()[i]=Data()[i]/a;
    return ret;
  }

  inline void operator/=(const MatrixBase<ROW,COL, T> &a) {
    for (unsigned int i=0; i<ROW*COL; i++)
      Data()[i]/=a.Data()[i];
  }

  inline void operator/=(const T &a) {
    for (unsigned int i=0; i<ROW*COL; i++)
      Data()[i]/=a;
  }

  // matrix math
  inline void Diagonal(const T x) {
    ZCHECK_EQ(ROW, COL);
    Zero();
    for (int i=0; i<ROW; i++)
      At(i, i)=x;
  }

  inline void Diagonal(const Vector<ROW, T> &x) {
    ZCHECK_EQ(ROW, COL);
    Zero();
    for (int i=0; i<ROW; i++)
      At(i, i)=x[i];
  }

  inline static MatrixBase<ROW, COL, T> GetDiagonal(const T x) {
    MatrixBase<ROW, COL, T> mat;
    mat.Diagonal(x);
    return mat;
  }

  inline static MatrixBase<ROW, COL, T> GetDiagonal(const Vector<ROW, T> &x) {
    MatrixBase<ROW, COL, T> mat;
    mat.Diagonal(x);
    return mat;
  }

  inline void Identical() {
    Diagonal(1);
  }

  inline static MatrixBase<ROW, COL, T> GetIdentical() {
    MatrixBase<ROW, COL, T> mat;
    mat.Identical();
    return mat;
  }

  inline void Transpose() {
    ZCHECK_EQ(ROW, COL);
    T t;
    for (int i=0; i<ROW; i++) for (int j=i; j<ROW; j++)
      if (i==j) continue;
      else {
        t=At(i, j);
        At(i, j)=At(j, i);
        At(j, i)=t;
      }
  }

  inline MatrixBase<COL,ROW, T> Transposed() const {
    MatrixBase<COL,ROW, T> ret;
    for (int i=0; i<ROW; i++) for (int j=0; j<COL; j++)
      ret.At(j, i)=At(i, j);
    return ret;
  }

  T Determinant() const {
    ZCHECK_EQ(ROW, COL);
    MatrixBase<ROW,COL, T> a = *this;
    int indx[ROW];  
    T d;
    a.luDecomposition(&indx[0], d);
    for(int i=0; i<ROW; i++) d *= a.V[i][i];      
    return d;
  }

  bool Invert() {
    ZCHECK_EQ(ROW, COL);
    MatrixBase<ROW,COL, T> a = *this;
    int indx[ROW];

    Vector<ROW, T> col;
    T d;
    int i, j;
    if (!a.luDecomposition(&indx[0], d)) //Decompose the matrix just once.
      return false; //singular matrix
    for(j=0; j<ROW; j++) { //Find inverse by columns.
      for(i=0; i<ROW; i++) col[i]=0.0;
      col[j]=1.0;
      a.luBackSubstitution(&indx[0],col);
      for(i=0; i<ROW; i++) 
        V[i][j]=col[i];    
    }
    return true;
  }

  MatrixBase<ROW,COL, T> Inverted() const {
    ZCHECK_EQ(ROW, COL);
    MatrixBase<ROW,COL, T> a = *this;
    MatrixBase<ROW,COL, T> y;
    int indx[ROW];

    Vector<ROW, T> col;
    T d;
    int i, j;
    a.luDecomposition(&indx[0], d); //Decompose the matrix just once.
    for(j=0; j<ROW; j++) { //Find inverse by columns.
      for(i=0; i<ROW; i++) 
        col[i]=0.0;
      col[j]=1.0;
      a.luBackSubstitution(&indx[0],col);
      for(i=0; i<ROW; i++) 
        y.V[i][j]=col[i];    
    }
    return y;
  }

  bool IsSymmetric() const {
    if (ROW!=COL) return false;
    for (int i=0; i<ROW; ++i) {
      for (int j=0; j<ROW; ++j) {
        if (V[i][j] != V[j][i]) return false;          
      }
    }
    return true;
  }

  T Trace() const {
    ZCHECK_EQ(ROW, COL);
    T t=0;
    for (int i=0; i<ROW; i++) t+=V[i][i];
    return t;
  }

  // 
  // Code borrowed from the Numerical Recipes in C++
  // 
  void luBackSubstitution(int indx[ROW], Vector<ROW, T>& b) {
    ZCHECK_EQ(ROW, COL);
    int i,ii=0,ip, j;
    T sum;

    for (i=0; i<ROW; i++) {
      ip=indx[i];
      sum=b[ip];
      b[ip]=b[i];
      if (ii != 0)
        for (j=ii-1; j<i; j++) sum -= V[i][j]*b[j];
      else if (sum != 0.0)
        ii=i+1;
      b[i]=sum;
    }
    for (i=ROW-1; i>=0; i--) {
      sum=b[i];
      for (j=i+1; j<ROW; j++) sum -= V[i][j]*b[j];
      b[i]=sum/V[i][i];
    }
  }

  bool luDecomposition(int indx[ROW], T& d) {      
    ZCHECK_EQ(ROW, COL);
    int i,imax, j, k;
    T big,dum,sum,temp;    

    Vector<ROW, T> vv;
    d=1.0;
    for (i=0; i<ROW; i++) {
      big=0.0;
      for (j=0; j<ROW; j++)
        if ((temp=fabs(V[i][j])) > big) big=temp;
      if (big == 0.0) return false; //singular matrix
      vv[i]=1.0/big;
    }
    for (j=0; j<ROW; j++) {
      for (i=0; i<j; i++) {
        sum=V[i][j];
        for (k=0; k<i; k++) sum -= V[i][k]*V[k][j];
        V[i][j]=sum;
      }
      big=0.0;
      for (i=j; i<ROW; i++) {
        sum=V[i][j];
        for (k=0; k<j; k++) sum -= V[i][k]*V[k][j];
        V[i][j]=sum;
        if ((dum=vv[i]*fabs(sum)) >= big) {
          big=dum;
          imax=i;
        }
      }
      if (j != imax) {
        for (k=0; k<ROW; k++) {
          dum=V[imax][k];
          V[imax][k]=V[j][k];
          V[j][k]=dum;
        }
        d = -d;
        vv[imax]=vv[j];
      }
      indx[j]=imax;
      if (V[j][j] == 0.0) V[j][j]=EPS; //_TINY;
      if (j != ROW-1) {
        dum=1.0/(V[j][j]);
        for (i=j+1; i<ROW; i++) V[i][j] *= dum;
      }
    }
    return true;
  }

    //c'tor and d'tor
  MatrixBase(){}
  explicit MatrixBase(const T &a){for (size_type i=0; i<ROW*COL; i++) Data()[i]=a;}
  template<typename T1> explicit MatrixBase(const T1 *data){for (size_type i=0; i<ROW*COL; i++) Data()[i]=data[i];}
  template<typename T1> explicit MatrixBase(const MatrixBase<ROW,COL,T1> &data){for (zuint i=0; i<ROW*COL; i++) Data()[i]=T(data[i]);}
  MatrixBase(const MatrixBase<ROW,COL, T> &other){*this=other;}
};

//multiply vector
template <zuint ROW,zuint COL,typename T>
inline const Vector<ROW, T> operator*(const MatrixBase<ROW,COL, T> &mat, const VectorBase<COL, T> &vec) {
  Vector<ROW, T> ret;
  for (Vector<ROW, T>::size_type i=0; i<ROW; i++)
  {
    ret[i]=T(0);
    for (Vector<ROW, T>::size_type j=0; j<COL; j++)
      ret[i]+=mat(i, j)*vec[j];
  }
  return ret;
}

template <zuint ROW,zuint COL,typename T>
inline const Vector<COL, T> operator*(const VectorBase<ROW, T> &vec,const MatrixBase<ROW,COL, T> &mat) {
  Vector<COL, T> ret;
  for (Vector<COL, T>::size_type i=0; i<COL; i++)
  {
    ret[i]=T(0);
    for (Vector<COL, T>::size_type j=0; j<ROW; j++)
      ret[i]+=vec[j]*mat(j, i);
  }
  return ret;
}

//scale at front
template <zuint ROW,zuint COL,typename T>
inline const MatrixBase<ROW,COL, T> operator+(const T &v,const MatrixBase<ROW,COL, T> &mat) {
  return mat+v;
}

template <zuint ROW,zuint COL,typename T>
inline const MatrixBase<ROW,COL, T> operator-(const T &v,const MatrixBase<ROW,COL, T> &mat) {
  return -mat+v;
}

template <zuint ROW,zuint COL,typename T>
inline const MatrixBase<ROW,COL, T> operator*(const T &v,const MatrixBase<ROW,COL, T> &mat) {
  return mat*v;
}

template <zuint ROW,zuint COL,typename T>
inline const MatrixBase<ROW,COL, T> operator/(const T &v,const MatrixBase<ROW,COL, T> &mat) {
  MatrixBase<ROW,COL, T> res;
  for (size_type i=0; i<ROW*COL; i++) res[i]=v/mat[i];
  return res;
}

//write all these stuff just to avoid direct memory copy when construct and copy
template<zuint ROW, zuint COL, typename T>
class Matrix : public MatrixBase<ROW,COL, T>
{
public:
  //constructor
  Matrix(void){}
  Matrix(const Matrix<ROW,COL, T>& a):MatrixBase<ROW,COL, T>(a) {}
  Matrix(const MatrixBase<ROW,COL, T>& a):MatrixBase<ROW,COL, T>(a) {}
  template<typename T1> explicit Matrix(const MatrixBase<ROW,COL,T1>& a):MatrixBase<ROW,COL, T>(a) {}
  template<typename T1> explicit Matrix(const T1 *p):MatrixBase<ROW,COL, T>(p) {}
  explicit Matrix(const T &a):MatrixBase<ROW,COL, T>(a){}
  template<typename T1> explicit Matrix(const Matrix<ROW,COL,T1>& a):MatrixBase<ROW,COL, T>(a) {}

  //assign
  inline const Matrix<ROW,COL, T> &operator=(const Matrix<ROW,COL, T>& a){MatrixBase<ROW,COL, T>::operator =(a); return *this;}
  template <typename T1>
  inline const Matrix<ROW,COL, T> &operator=(const Matrix<ROW,COL,T1>& a){MatrixBase<ROW,COL, T>::operator =(a); return *this;}
};

// IOObject
template<unsigned int ROW, unsigned int COL, typename T>
class IOObject<MatrixBase<ROW,COL, T> >
{
public:
  static void WriteFileB(FILE *fp, const MatrixBase<ROW,COL, T>* src, zsize size) {
    IOObject<T>::WriteFileB(fp, src->Data(), size * ROW * COL);
  }
  static void ReadFileB(FILE *fp, MatrixBase<ROW,COL, T>* dst, zsize size) {
    IOObject<T>::ReadFileB(fp, dst->Data(), size * ROW * COL);
  }
  static void WriteFileR(RecordFile &fp, const zint32 label, const MatrixBase<ROW,COL, T>* src, zsize size) {
    IOObject<T>::WriteFileR(fp, label, src->Data(), size * ROW * COL);
  }
  static void ReadFileR(RecordFile &fp, const zint32 label, MatrixBase<ROW,COL, T>* dst, zsize size) {
    IOObject<T>::ReadFileR(fp, label, dst->Data(), size * ROW * COL);
  }
  static void WriteFileR(RecordFile &fp, const zint32 label, const MatrixBase<ROW,COL, T>& src) {
    IOObject<MatrixBase<ROW,COL, T> >::WriteFileR(fp, label, &src, 1);
  }
  static void ReadFileR(RecordFile &fp, const zint32 label, MatrixBase<ROW,COL, T>& dst) {
    IOObject<MatrixBase<ROW,COL, T> >::ReadFileR(fp, label, &dst, 1);
  }
  static void CopyData(MatrixBase<ROW,COL, T>* dst, const MatrixBase<ROW,COL, T>* src, zsize size) {
    IOObject<T>::CopyData(static_cast<T*>(dst->Data()), static_cast<const T*>(src->Data()), size*ROW*COL);
  }
};

template<unsigned int ROW, unsigned int COL, typename T>
class IOObject<Matrix<ROW,COL, T> >
{
public:
  static void WriteFileB(FILE *fp, const Matrix<ROW,COL, T>* src, zsize size) {
    IOObject<MatrixBase<ROW,COL, T> >::WriteFileB(fp, src, size);
  }
  static void ReadFileB(FILE *fp, Matrix<ROW,COL, T>* dst, zsize size) {
    IOObject<MatrixBase<ROW,COL, T> >::ReadFileB(fp, dst, size);
  }
  static void WriteFileR(RecordFile &fp, const zint32 label, const Matrix<ROW,COL, T>* src, zsize size) {
    IOObject<MatrixBase<ROW,COL, T> >::WriteFileR(fp, label, src, size);
  }
  static void ReadFileR(RecordFile &fp, const zint32 label, Matrix<ROW,COL, T>* dst, zsize size) {
    IOObject<MatrixBase<ROW,COL, T> >::ReadFileR(fp, label, src, size);
  }
  static void WriteFileR(RecordFile &fp, const zint32 label, const Matrix<ROW,COL, T>& src) {
    IOObject<MatrixBase<ROW,COL, T> >::WriteFileR(fp, label, src);
  }
  static void ReadFileR(RecordFile &fp, const zint32 label, Matrix<ROW,COL, T>& dst) {
    IOObject<MatrixBase<ROW,COL, T> >::ReadFileR(fp, label, dst);
  }
  static void CopyData(Matrix<ROW,COL, T>* dst, const Matrix<ROW,COL, T>* src, zsize size) {
    IOObject<MatrixBase<ROW,COL, T> >::CopyData(dst, src, size);
  }
};


template <zuint ROW, zuint COL, typename T> inline Matrix<ROW, COL, T> Abs(const MatrixBase<ROW, COL, T> &x) {
  Matrix<ROW, COL, N> ret;
  for (size_type i=0; i<ret.size(); i++) ret[i]=Abs(x[i]);
  return ret;
}


}
